%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code is written to estimate the GARCH-MIDAS model with fixed span RV
%
% This code needs the data of the following format:
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stkmkt_ret_data
% yyyy  mm  dd  yyyymmdd  vwretd  ewretd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% For practical purposes, the parameter m and theta in the tau function
% are squared as follows:
%
% tau = (m^2) + (theta^2) * ( X * betapolyn(K,[(K-1):-1:1]',w) )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

% Determine [fix] and [K]--------------------------------------------------------
fix  = input('How long should tau be fixed? (monthly/quarterly/biannually/annually)   ','s');
disp(' ')
K    = input('How many lags do you want to include in MIDAS?     ');
% -------------------------------------------------------------------------------

% Determine [fix_n] accordingly -------
switch fix
    case 'monthly'
        fix_n=1;
    case 'quarterly'
        fix_n=3;
    case 'biannually'
        fix_n=6;
    case 'annually'
        fix_n=12;
    otherwise
        disp('wrong answer')
        stop
end;
% ------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% daily stock returns: 1926 - 2007 / NYSE-AMEX-NASDAQ
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stkmkt_ret_data
% yyyy  mm  dd  yyyymmdd  vwretd  ewretd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load stkmkt_ret_data;
daily_stk = [stkmkt_ret_data(:,1:3) stkmkt_ret_data(:,5)];


[T0,C0]=size(daily_stk);

ind0=find( diff(daily_stk(:,2)) ~= 0 ); % row number
ind_mth=zeros(T0,1);
ind0=ind0+1; % row number where new month begins
             % ex) [15;45;75;105;135;165;...]'

%-------------
ind_mth(ind0,1)=1; 
%-------------

daily_stk=[daily_stk(:,1:3) ind_mth daily_stk(:,4)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine sample period 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

date=daily_stk(:,1)*10000 + daily_stk(:,2)*100 + daily_stk(:,3);

disp('    ')
disp('Sample Period? (ex. [19270101;20071231] / for whole sample, just press enter)');
disp('The first (fixed period) X (number of lags) data will be burned');
subperiod=input('Your Choice?   ');

if isempty(subperiod)==1
    subperiod=[date(1);date(end)];
else
end;

start_ind=find(date==subperiod(1));
j=1;
while isempty(start_ind)==1
    start_ind=find(date==(subperiod(1)+j));
    j=j+1;
end;

end_ind=find(date==subperiod(2));
j=1;
while isempty(end_ind)==1
    end_ind=find(date==(subperiod(2)-j));
    j=j+1;
end;

daily_stk=daily_stk(start_ind:end_ind,:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adding timely index and number of days in the period %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T0,C0]=size(daily_stk);

ft_ind0=find(daily_stk(:,4)==1); % [1;31;61;91;121;151;181;...]'
L=length(ft_ind0);  % number of months included in the sample
fix_v=[1:fix_n:L]'; % [1;4;7;10;13;....]'
ft_ind0=ft_ind0(fix_v);  % fixed term month change index (0 or 1) to be included in daily_stk data
                         % [1;91;181;.....]' corresponding to {4/1,7/1,...)
L2=length(ft_ind0);
ft_ind=zeros(T0,1);
ft_ind(ft_ind0)=1; % assigning 1 to 4/1,7/1,...

%--------------------------------------------------------------
daily_stk=[daily_stk(:,1:4) ft_ind daily_stk(:,5)];
%--------------------------------------------------------------
ind1=find( diff(daily_stk(:,2)) ~= 0 ); % row number
ind1=ind1+1; %

num_days = zeros(T0,1);
num_days(ind1)=[diff(ind1);T0-ind1(end)];
num_days(1) = ind1(1)-1;
%--------------------------------------------------------------
daily_stk=[daily_stk(:,1:5) num_days daily_stk(:,6)];
%--------------------------------------------------------------

num_days_fx      =zeros(T0,1);
num_days_fx_temp =daily_stk(:,4).*daily_stk(:,6);

for i=1:(L2-1)
    num_days_fx(ft_ind0(i))=sum( num_days_fx_temp( ft_ind0(i):(ft_ind0(i+1)-1) ) );
end;

num_days_fx( ft_ind0(L2) )=sum( num_days_fx_temp( ft_ind0(L2):end ) );

%--------------------------------------------------------------
daily_stk=[daily_stk(:,1:6) num_days_fx daily_stk(:,7)];
%--------------------------------------------------------------
% (1) yr 
% (2) mth 
% (3) day 
% (4) 'beginning of new month index (0 or 1)' 
% (5) 'beginning of new fixed term (0 or 1)'
% (6) num_days
% (7) num_days_fx
% (8) ret
%-------------------------
ave_days_fx=mean(daily_stk( find(daily_stk(:,7)~=0),7 ));

r=daily_stk(:,8);
r2=r.^2;

%-------------------------
N_ft=sum(daily_stk(:,5));
RV=zeros(N_ft,1);
skew=zeros(N_ft,1);

%--------------------------------------------------
for i=1:(N_ft-1)
    RV(i)    = sum( r2( ft_ind0(i):(ft_ind0(i+1)-1) ) ); % / daily_stk(ft_ind0(i),7) * ave_days_fx;
end;

RV(N_ft)   = sum( r2( ft_ind0(N_ft):end ) ); % / daily_stk(ft_ind0(end),7) * ave_days_fx;
        
           
% RV is realized volatility over the fixed term
%-------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setting new "daily_stk" starting from the next of lagged terms 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
daily_stk=daily_stk(ft_ind0(K+1):end,:);  %(1) changed K -> 48/16/8
r=daily_stk(:,8);
r2=r.^2;

start_yr  = daily_stk(1,1);
end_yr    = daily_stk(end,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

T0=length(r);

ft_ind_from=zeros(T0,1);
ft_ind_from(1)=K-K+1; %(2) changed 1 -> 10-K+1
ft_ind_to=zeros(T0,1);
ft_ind_to(1)=K;   %(3) changed K -> 10

for k=2:T0
    ft_ind_from(k) = ft_ind_from(k-1) + daily_stk(k,5);  % start row # of fixed term average X
    ft_ind_to(k)   = ft_ind_to(k-1)   + daily_stk(k,5);      % end row # of fixed term average X
end;

ft_ind_from=ft_ind_from + 1; %recall that we're taking only K-1 lagged terms

%--------------------------------------------------------------
daily_stk=[daily_stk(:,1:7) ft_ind_from ft_ind_to daily_stk(:,8)];
% (1) yr 
% (2) mth 
% (3) day 
% (4) 'beginning of new month index (0 or 1)' 
% (5) 'beginning of new fixed term (0 or 1)'
% (6) num_days
% (7) num_days_fx
% (8) 'from-row-index' to take data from 'monthly_rate'
% (9) 'to-row-index' to take data from 'monthly_rate'
% (10) ret
%-----------------------------------------------------

disp('starting  date for r')
daily_stk(1,1:3)

% Constructing Quarterly RV series -------------------------
qtr_ind = find( daily_stk(:,5)==1 );
if qtr_ind(1)~=1
    qtr_ind = [1;qtr_ind];
else
end;

N_qtr  = length(qtr_ind);
RV_qtr = zeros(N_qtr,1);
Quarticity = zeros(N_qtr,1);

for i=1:(N_qtr-1)
    RV_qtr(i) = sum( daily_stk( qtr_ind(i):(qtr_ind(i+1)-1) , 10 ).^2 );
    Quarticity(i) = sum( daily_stk( qtr_ind(i):(qtr_ind(i+1)-1) , 10 ).^4 );
end;

RV_qtr(end) = sum( daily_stk( qtr_ind(N_qtr):end , 10 ).^2 );
%--------------------------------------------------------



% constructing daily MIDAS-RV series ------------------
X=zeros(T0,K-1);
for j=1:T0
    ft_ind_from=daily_stk(j,8);
    ft_ind_to=daily_stk(j,9);
    X(j,:)=RV(ft_ind_from:ft_ind_to)';
end;
%-----------------------------------------------------

RV=RV(K+1:end); %(4)

r_yr = daily_stk(:,1);

ref_yr = [start_yr+(5-mod(start_yr,5)):5:end_yr - mod(end_yr,5)]';

if mod(start_yr,5)==0
    ref_yr=[start_yr;ref_yr];
else
end;
        
ref_yr_length=length(ref_yr);
tick_yr=zeros(ref_yr_length,1);

for i=1:ref_yr_length
    tick_yr(i)=min(find(r_yr==ref_yr(i)));
end;



%-----------------------------------------------------
% ESTIMATION
%-----------------------------------------------------

miu_i    = 0.0005;
alpha_i  = 0.08;
beta_i   = 0.90;
theta_i  = 0.1;
w_i      = 5;
m_i      = 0.0001;

parameter_i=[miu_i;alpha_i;beta_i;theta_i;w_i;m_i];

g1=1;

options=optimset('Display','iter','MaxFunEvals',3000,'MaxIter',3000,'TolFun',1e-6,'TolX',1e-6);
options =optimset(options,'GradObj','off');
options=optimset(options,'DerivativeCheck','off');
options=optimset(options,'Diagnostics','off');

A=[
    0  1  1  0  0  0;
    0 -1  0  0  0  0;
    0  0 -1  0  0  0;
    0  0  0  0 -1  0;
    0  0  0  0  1  0;
    0  0  0  0  0 -1;
  ];

b=[
 .999999999;
          0; 
          0;
         -1;
       1000;
 -.00000001
   ];

   
disp('Optimization Results ------------------')
[parameter,fval,exitflag,output,lambda,grad,hessian] = fmincon('nlh_garchmidas',parameter_i,A,b,[],[],[],[],[],options,g1,K,r,X)
disp('---------------------------------------')

T_n = length(r);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cov_Var = hessian/T_n;
V       = inv(Cov_Var);
S_var   = diag(V)/T_n;
S_err   = sqrt(S_var);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

LLF=-nlh_garchmidas(parameter,g1,K,r,X);
disp(['Likelihood= ' num2str(LLF)])
disp(' ')

num_par = 6;

% BIC = -2*LLF/T_n + num_par*log(T_n)/T_n

name=['miu=   ';'alpha= ';'beta=  ';'theta= ';'w=     ';'m=     '];

disp('name      estimates        ')
disp([name num2str(parameter)])
disp('   ')
disp('name      s.e.        ')
disp([name num2str(S_err)])
disp('   ')
disp('name      t-stat        ')
tstat=parameter./S_err;
disp([name num2str(tstat)])
disp(' ')

miu   =parameter(1);
alpha =parameter(2);
beta  =parameter(3);
theta =parameter(4);
w     =parameter(5);
m     =parameter(6);

%---------------------------------------------------------------
tau = (m^2) + (theta^2) * ( X * betapolyn(K,[(K-1):-1:1]',w) ); % function of theta and w 
%---------------------------------------------------------------

%--------------------------------------------------------------------
g=zeros(T_n,1);
g(1)=g1;

for i=2:T_n
    g(i)=(1-alpha-beta)+alpha*((r(i-1)-miu)^2)/tau(i-1)+beta*g(i-1);
end;
%--------------------------------------------------------------------

e=(r-miu)./sqrt(tau.*g);

figure
hist(e,100);title('Distribution of errors')

figure
plot(flipdim(betapolyn(K,[(K-1):-1:1]',w),1));title('optimal MIDAS weights')



figure
plot(sqrt(252*tau.*g),'g--','LineWidth',1);
hold on
plot(sqrt(252*tau),'b-','LineWidth',2);
axis tight
legend('(ann.) conditional volatility (\tau*g)^{1/2}','(ann.) secular component (\tau)^{1/2}')
title('conditional volatility and its long run component of stock market returns');
xlabel('year')
ylabel('annualized volatility')
set(gca,'XTick',tick_yr);
set(gca,'XTickLabel',num2str(ref_yr));


